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The isotropic to nematic transition in a system of soft spherocylinders is studied by means of grand 
canonical Monte Carlo simulations. The probability distribution of the particle density is used to 
determine the coexistence density of the isotropic and the nematic phases. The distributions are 
also used to compute the interfacial tension of the isotropic-nematic interface, including an analysis 
of finite size effects. Our results confirm that the Onsager limit is not recovered until for very large 
elongation, exceeding at least L/D = 40, with L the spherocylinder length and D the diameter. For 
smaller elongation, we find that the interfacial tension increases with increasing L/D, in agreement 
with theoretical predictions. 
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I. INTRODUCTION 

On change of density, suspensions of rod-like parti- 
cles undergo a phase transition between an isotropic fluid 
phase, where the particle orientations are evenly dis- 
tributed, and an anisotropic nematic fluid phase, where 
the particle orientations are on average aligned. This 
phenomenon was explained by Onsager in a theory based 
on infinitely elongated hard spherocylinders Q . Onsager 
theory has been remarkably successful at describing the 
isotropic to nematic (IN) transition, and still serves as 
the basis for many theoretical investigations of the prop- 
erties of liquid crystals. Over the last twenty years, for 
instance, several groups have investigated the properties 
of the IN interface usingOnsager-type density functional 
approaches 0i0,SSIail3 ■ An important finding of these 
studies is that the interfacial tension tin of the IN inter- 
face is minimized when the director, which is the axis 
of average orientation of the particles, lies in the plane 
of the interface. In the case of in-plane alignment, tin 
is predicted to be very low, but the precise value varies 
considerably between different authors Theoretical 
estimates for tin typically range from 0.156 Q to 0.34 
0, in units of k B T/LD, with L the rod length, D the 
rod diameter, T the temperature, and k B the Boltzmann 
constant. 

Obviously, the Onsager limit of infinite rod length is 
purely academic. In order to describe more realistic situ- 
ations, it is necessary to go beyond the Onsager approx- 
imation, and consider the case of finite rod length. An 
example is the theoretical work of Ref. El which demon- 
strates that the interfacial tension in the case of finite rod 
length is considerably lower than predicted by Onsager 
theory. 

To test the accuracy of the theoretical estimates of tin , 
one might envision a direct comparison to experimental 
data. Unfortunately, this is not straightforward. The 
models used in theoretical treatments of the IN interface 
are typically rather simplistic, usually based on a short- 
ranged pair potential in a system of monodisperse sphe- 
rocylinders. It is not reasonable to expect quantitative 



agreement with experiments using these models, because 
the interactions in the experimental system will be much 
more complex. For example, polydispersity may be an 
important factor, and it is not clear to what extent long- 
range interactions play a role. Even the experimental 
determination of the rod dimensions L and D, required 
if a comparison to theory is to be made, presents compli- 
cations |j- 

In order to validate the assumptions made by the 
various theoretical approaches, it is nevertheless im- 
portant to test the accuracy of the theoretical predic- 
tions. To this end, computer simulations are ideal, 
because they, in principle, probe the phase behavior 
of the model system without resorting to approxima- 
tions. With inexpensive computer power readily avail- 
able nowadays, several groups have taken the opportu- 
nity to investigate the IN transition by means of simula- 
tions [H E El H El 13 IH El An example of this 
approach is Ref. [lj, where the coexistence properties of 
the bulk isotropic and nematic phases of hard sphero- 
cylinders are carefully mapped out using Gibbs ensemble 
Monte Carlo El- These simulations generally recover 
the Onsager limit for long rods, while for shorter rods 
pronounced deviations show up jl2j |. Unfortunately, the 
Gibbs ensemble cannot be used to measure tin 7 which is 
the aim of this work. 

To obtain tin in simulations, different techniques must 
be used. One such technique is based on the anisotropy of 
the pressure tensor. In Ref. ll8l this method is applied to 
suspensions of ellipsoids with axial ratio k = A/ B = 15, 
where A is the length of the symmetry axis of the el- 
lipsoids, and B that of the transverse axis. The corre- 
sponding interfacial tension is 0.006 ± 0.005 k B T/B 2 ss 
0.09 k B T/AB if a hard interaction potential is used, and 
0.011 ± 0.004 k B T/B 2 w 0.165 k B T/AB using a soft po- 
tential. Note that the anisotropy of the pressure tensor is 
very small, and therefore difficult to measure accurately 
in practice, as indicated by the error bars. 

In Ref. 123, again for (soft) ellipsoids with k = 
15, a value of the interfacial tension tin = 0.016 ± 
0.002 k B T/B 2 fa 0.24 k B T/AB is reported. This re- 
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suit was obtained by measuring the capillary broadening 
of the IN interface. According to capillary wave theory 
[2ll | , the mean squared amplitudes of the capillary fluctu- 
ations are proportional to I/71N, and this can be used to 
obtain the interfacial tension. Unfortunately, capillary 
wave theory is only valid in the long wavelength limit, 
such that very large system sizes are required. Moreover, 
if, as in Ref. l20l periodic boundary conditions are used, 
two interfaces will be present in the simulation box. Since 
7in is very small, large capillary fluctuations can occur, 
and one needs to be aware of interactions between the 
two interfaces. 

Clearly, in order to obtain 71N more accurately, much 
more computer power or different simulation techniques 
are required. Recent advances in grand canonical sam- 
pling methods [22I |23| have enabled accurate measure- 
ments of the interfacial tension in simple fluids E3, |2j| , 
and complex fluids such as polymer solutions |26( and 
colloid-polymer mixtures |27j . The aim of this paper is 
to apply these techniques to the IN transition in a system 
of soft spherocylinders, and to extract the corresponding 
phase diagram and the interfacial tension. Simulations 
in the grand canonical ensemble offer a number of ad- 
vantages over the more conventional methods discussed 
previously. More precisely, in grand canonical simula- 
tions, both the coexistence properties can be probed, as 
in the Gibbs ensemble, as well as the interfacial proper- 
ties. Additionally, finite - size scaling methods are avail- 
able which can be used to extrapolate simulation data 
to the thermodynamic limit |28l 123 . l3Ct |3l| . It has been 
demonstrated that grand canonical ensemble simulations 
combined with novel finite size scaling algorithms can 
yield results of truly impressive accuracy 31,]. 

This article is structured as follows: First, we introduce 
the soft spherocylinder model used in this work. Next, 
we describe the grand canonical Monte Carlo method, 
and explain how the coexistence properties, and the in- 
terfacial tension are obtained. Finally, we present our 
results, followed in the last section by a discussion and 
an outlook to future work. 



II. MODEL 

In this study, the particles are modeled as repulsive 
soft spherocylinders of elongation L and diameter D. For 
numerical convenience, a very simple potential has been 
chosen. The interaction between two rods A and B, is 
given by a pair potential of the form 



Vab(t) 



e r < D, 
otherwise, 



(1) 



with r the distance between two line segments of length 
L, see Fig. 2] The total energy is thus proportional to the 
number of overlaps in the system. In this work, the rod 
diameter D is taken as unit of length, and k^,T as unit 
of energy. The strength of the potential is set to e = 2. 
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FIG. 1: Two-dimensional representation of the simulation 
model of this work. The liquid crystals are modeled as soft 
spherocylinders with elongation L and diameter D. Two rods 
A and B interact via the pair potential of Eq.Q, which is a 
function of their minimum distance r only. If the rods overlap, 
the system pays a constant energy cost e. To speed up the 
determination of overlap, the simulation box is subdivided 
into cubic cells with edge length a, see details in text. 



Note that in the limit e — > 00, this model approaches a 
system of infinitely hard rods. 

To study the IN transition, we typically use the den- 
sity and the rod alignment as order parameters. Note 
that both the isotropic and the nematic phase are fluid 
phases, in the sense that long-range positional order of 
the centers of mass is absent. In the nematic phase, how- 
ever, there is orientational order where, on average, the 
rods point in one direction (called the director). In the 
isotropic phase, on the other hand, there is no orienta- 
tional order. Since the density of the nematic phase is 
slightly higher than of the isotropic phase, we may use 
the particle number density p = N/V to distinguish be- 
tween both phases, with N the number of rods in the 
system, and V the volume of the simulation box. Fol- 
lowing convention, we also introduce the reduced density 
P* = p/Pc P , with p cp = 2/(V2+ (L/D)VS) the density 
of regular close packing of hard spherocylinders. Ori- 
entational order is as usual measured by the S2 order 
parameter, defined as the maximum eigenvalue of the 
orientational tensor Q: 



1 N 



(2) 



Here, Ui a is the a component (a = x, y, z) of the orien- 
tation vector Ui of rod i (normalized to unity), and S a p 
is the Kronecker delta. In case of orientational order, 
such as in the nematic phase, 5 2 assumes a value close to 
one, while in the disordered isotropic phase, S2 is close 
to zero. 
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FIG. 2: Coexistence distribution W = k B T In P(N) of the 
isotropic to nematic transition in a system of soft rods inter- 
acting via Eq.Q with e = 2 and L/D = 15. The low density 
peak corresponds to the isotropic phase (ISO), the high den- 
sity peak to the nematic phase (NEM) , and the barrier AF to 
the free energy difference between the two phases (AF is given 
by the average peak height as measured from the minimum 
in between the peaks). The above distribution was obtained 
using box dimensions L x = 2.1L and L z = SAL. The coexis- 
tence value of the chemical potential reads [i = 5.15 and was 
obtained using the equal area criterion described in the text. 



III. SIMULATION METHOD 

The simulations are performed in the grand canoni- 
cal ensemble. In this ensemble, the volume V, the tem- 
perature T, and the chemical potential /i of the rods 
are fixed, while the number of rods N inside the sim- 
ulation box fluctuates. Insertion and removal of rods 
are attempted with equal probability, and accepted with 
the standard grand canonical Metropolis rules, given by 



A(N N + 1) 



N-l) 



iin 1 . 

-/3AE- 



N+V 
/3 M 



and A(N 



] , with AE the energy dif- 



ference between initial and final state, and (3 — 1/ksT 
|29ll32j |. The simulations are performed in a three dimen- 
sional box of size L x x L y x L z using periodic boundary 
conditions in all directions. In this work, we fix L x = L y , 
but we allow for elongation in the remaining direction 
Lz > L x . Moreover, to avoid double interactions between 
rods through the periodic boundaries, we set L x > 2L. 

During the simulations, we measure the probability 
P(N), defined as the probability of observing a system 
containing TV rods. Note that the shape of the dis- 
tribution will depend on the rod elongation L/D, the 
temperature T, and the chemical potential fi. More- 
over, there may be finite-size effects, introducing addi- 
tional dependences on the box dimensions L x and L z . 
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FIG. 3: Snapshot of a system of soft spherocylinders at IN co- 
existence. The spherocylinders are shaded according to their 
orientation. On the left side of the dashed line the system is 
isotropic, on the right side it is nematic. The second inter- 
face coincides with the boundaries of the box in the elongated 
direction. 



At phase coexistence, the distribution P(N) becomes 
bimodal, with two peaks of equal area, one located at 
small values of N corresponding to the isotropic phase, 
and one located at high values of N corresponding to 
the nematic phase. A typical coexistence distribution is 
shown in Fig. where the logarithm of P(N) is plot- 
ted. Coexistence is determined using the equal area 
rule |33j | . At coexistence, the equal area rule implies 
that / <JV> P{N)dN = /~ P(N)dN, with (N) the av- 
erage of the full distribution (N) = NP(N)dN, 
where we assume that P(N) has been normalized to 
unity J Q P(N)dN = 1. The coexistence density of the 
isotropic phase follows trivially from the average of P(N) 

in first peak pi SO = (2/V) J^ N) NP(N)dN, and similarly 
for the nematic phase pnem = (2/V) J*,^ NP(N)dN, 
where the factors of two are a consequence of the nor- 
malization of P(N). 

The interfacial tension 7in is extracted from the loga- 
rithm of the probability distribution W = kBT\nP(N). 
Since — W corresponds to the free energy of the system, 
the average height AF of the peaks in W, measured with 
respect to the minimum in between the peaks, equals the 
free energy barrier separating the isotropic from the ne- 
matic phase. When the overall density of the system is 
in the interval between the peaks piso < p < Pnem, 
coexistence between an isotropic and nematic domain is 
observed. A snapshot of the system in this regime re- 
veals a slab geometry, with one isotropic region, and one 
nematic region, separated by an interface (because of pe- 
riodic boundary conditions, there are actually two inter- 
faces). An example snapshot is shown in Fig. Note 
that the director of the nematic phase lies in the plane 
of the interfaces. This was the typical case for the snap- 
shots studied by us, and is consistent with the theoretical 
prediction that in-plane alignment yields the lowest free 
energy. 

The barrier AF in Fig. thus corresponds to the free 
energy cost of having two interfaces in the system. Since, 
in this work, the box dimensions are chosen such that 
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L x = L y and L z > L x , the interfaces will be oriented 
perpendicular to the elongated direction, since this min- 
imizes the interfacial area, and hence the free energy of 
the system. The total interfacial area in the system thus 
equals 2L\. Since the interfacial tension is simply the 
excess free energy per unit area, we may write 

7 in(Lx) = AF/(2L 2 X ), (3) 

with Tin (Ac) the interfacial tension in a finite simula- 
tion box with lateral dimension L x |2q]. To obtain the 
interfacial tension in the thermodynamic limit, one can 
perform a finite size scaling analysis to estimate 
limi x _^ 00 7in(^x)- Alternatively, away from any criti- 
cal point, the most dominant finite size effects will likely 
stem from interactions between the two interfaces. In 
this case, it is feasible to use an elongated simulation 
box with L z 3> i x , such as in Fig. |3J The advantage 
of using an elongated simulation box is that interactions 
between the interfaces are suppressed. This enhances a 
flat region in W between the peaks, indicating that the 
interfaces are no longer interacting, and that finite size 
effects will likely be small. In this work, both approaches 
will be used. 

If the free energy barrier AF is large, transitions be- 
tween the isotropic and the nematic phase become less 
likely, and the simulation will spend most of the time 
in only one of the two phases. A crucial ingredient in 
our simulation is therefore the use of a biased sampling 
technique. We use successive umbrella sampling 23] to 
enable accurate sampling in regions where P(N), due 
to the free energy barrier separating the phases, is very 
small. Note also that phase coexistence is only observed 
if the chemical potential /i is set equal to its coexistence 
value. This value is in general not known at the start of 
the simulation, but it may easily be obtained by using the 
equation P(JV|/ii) = P(7V|p )e ,3(A11 ~ w,)Ar , with P(N\fj, a ) 
the probability distribution P(N) at chemical potential 
fi a . In the simulations, we typically set the chemical po- 
tential to zero and use successive umbrella sampling to 
obtain the corresponding probability distribution. We 
then use the above equation to obtain the desired coex- 
istence distribution, in which the area under both peaks 
is equal. 

IV. NUMERICAL OPTIMIZATIONS 

Most of the CPU time in our simulations is spend on 
calculating the distance r between two line segments, see 
Fig. ^ Naturally, one tries to minimize the number of 
calls to the routine that determines the distance. To 
this end, we use a cubic linked cell structure, which is 
schematically illustrated in Fig. ^ The crucial point is 
that the lattice constant a is chosen such that D < a < L. 
To determine if rod B in Fig. ^overlaps with any of the 
other rods in the system, it is sufficient to consider only 
those rods contained in the cubes intersected by rod B 
(shaded gray), plus the rods contained in the nearest and 



next-nearest neighbors of these cubes. Since the isotropic 
to nematic transition occurs at low density, most cubes 
will be empty, resulting in a substantial efficiency gain. 
Some CPU time is used for manipulating the linked cell 
structure, but for large systems (« N > 1500) and long 
rods (a* L/D > 10), the gain in efficiency is already a 
factor of five. Some fine-tuning is required to obtain 
an optimal value of the lattice constant. We found that 
a w 0.2L typically gives good results. 

A further optimization concerns the calculation of the 
S2 order parameter, see Eq.(|5J). In a naive implementa- 
tion, determining the orientational tensor Q involves an 
O(N) loop over all rods in the system. In our implemen- 
tation, the tensor elements of Q are updated after each 
accepted Monte Carlo move, which can be done at the 
cost of only a few additions and multiplications. Since we 
keep the tensor elements updated throughout the simu- 
lation, the 0{N) loop of Eq.(0| never needs to be carried 
out. Finally, to determine the maximum eigenvalue of Q, 
we do not use a numerical scheme, but instead use the 
exact expression for the roots of a third degree polyno- 
mial. The advantage of this implementation is that the 
value of 5*2 is known exactly throughout the simulation, 
at a cost exceeding no more than one percent of the total 
invested CPU time. 

We conclude this section with a few benchmarks. For 
e = 2 in Eq.QJ, we found that the acceptance rate 
of grand canonical insertion is around 9 percent in the 
isotropic phase, and it decreases to around 6 percent in 
the nematic phase. The acceptance rates are rather in- 
sensitive to L/D. With the optimized implementation 
described in this section, we can typically generate 5000- 
8000 accepted grand canonical moves per second on a 2.2 
GHz AMD Opteron processor. 



V. RESULTS 

A. Phase diagram 

We first use our grand canonical Monte Carlo scheme 
to determine the IN phase diagram of the soft sphero- 
cylinder system of Eq.QJ using e = 2. For several rod 
elongations L/D, we measured the distribution P(N), 
from which piso and pnem were obtained. The system 
size used in these simulations is typically L x = L y = 2.1L 
and L z = A.2L. In Fig.0] we plot the reduced density of 
the isotropic and the nematic phase as function of L/D. 
We observe that the phase diagram is qualitatively simi- 
lar to that of hard spherocylinders |l2( • The quantitative 
difference being that, for soft rods, the IN transition is 
shifted towards higher density. The inset of Fig. 0] shows 
the concentration variable c = irDL 2 p/4 as a function 
of D/L. For hard spherocylinders, Onsager theory pre- 
dicts that ciso = 3.29 and cnem = 4.19 in the limit of 
infinite rod length, or equivalently D/L — > 0. In case of 
the soft potential of Eq.QJ, these values must be mul- 
tiplied by (1 — e _/3e ) _1 w 1.16 for e — 2. In the inset 
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FIG. 4: Soft spherocy Under phase diagram of the IN transi- 
tion using e = 2 in Eq. Q . Shown is the reduced density p* of 
the isotropic phase (closed circles) and of the nematic phase 
(open circles) as function of L/D. The inset shows the con- 
centration variable c as function of D/L for both the isotropic 
and the nematic phase. The lower and upper arrow in the in- 
set mark the Onsager limit D/L — > for the isotropic and the 
nematic phase, respectively. The lines connecting the points 
serve as a guide to the eye. 



of Fig. ^ the corresponding limits are marked with ar- 
rows. As in Ref. 0, we observe that the simulation data 
for the isotropic phase smoothly approach the Onsager 
limit, while the nematic branch of the binodal seems to 
overshoot the Onsager limit. This we attribute to equi- 
libration problems. To simulate the IN transition in the 
limit D/L — > 0, large system sizes are required, and it be- 
comes increasingly difficult to obtain accurate results. To 
quantify the uncertainty in our measurements, additional 
independent simulations for rod elongation L/D = 25, 
30, and 35 were performed. The corresponding data are 
also shown in Fig. For L/D > 30, we observe sig- 
nificant scatter, while for L/D < 25, the uncertainty is 
typically smaller than the symbol size used in the plots. 



B. Interfacial tension 

Next, the interfacial tension just of the IN interface is 
determined for L/D = 10 and L/D — 15. Unfortunately, 
the system size used to compute the phase diagram in 
the previous section, was insufficient to accurately ex- 
tract the interfacial tension because no flat region be- 
tween the peaks in P(N) could be distinguished. This 
indicates that the interfaces are still strongly interacting. 
To properly extract the interfacial tension, much larger 
systems turned out to be required. In this case, care 
must be taken in the sampling procedure. Many sam- 
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FIG. 5: Monte Carlo time series of a biased grand canonical 
simulation. The top frame shows the S2 order parameter as 
a function of the invested CPU time, the lower frame the 
reduced density, with CPU time expressed in hours on a 2.6 
GHz Pentium. During the simulation, the reduced density 
was confined to the interval 0.245 < p* < 0.275, as indicated 
by the horizontal lines in the lower figure. The data were 
obtained using L/D = 15, e = 2, L x = 2.1L and L z = 8.4L, 
which are the same parameters as used in Fig. 



pling schemes, especially the ones that are easy to im- 
plement such as successive umbrella sampling, put a bias 
on the density only. Such schemes tend to "get stuck" 
in meta-stable droplet states when the system size be- 
comes large |2(j. As a result, one may have difficulty 
reaching the state with two parallel interfaces, in which 
case Eq.© cannot be used. 

Therefore, for large systems, one must carefully check 
the validity of the simulation results. To this end, we oc- 
casionally inspect simulation snapshots. For sufficiently 
elongated simulation boxes L z 3> L x and at densities in- 
side the coexistence region pi§o <C p <§C Pnem, we indeed 
observe two planar interfaces oriented perpendicular to 
the elongated direction, in accord with Fig. |3] To fur- 
ther check the consistency of the measured distributions 
P(N), we performed a number of additional grand canon- 
ical simulations using a biased Hamiltonian of the form 
H = Ho + W, with Ho the Hamiltonian of the real sys- 
tem defined by Eq.© and W = -k B T\nP(N). If the 
measured P(N) is indeed the equilibrium coexistence dis- 
tribution of the real system, a simulation using the biased 
Hamiltonian should visit the isotropic and the nematic 
phase equally often on average This is illustrated 

in the top frame of Fig. [SI which shows the S2 order 
parameter as a function of the elapsed simulation time 
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FIG. 6: Coexistence distributions W = k B T In P(N) of soft 
spherocylinders with L/D = 10 and e = 2 for various sys- 
tem sizes. In each of the above distributions, the lateral box 
dimension was fixed at L x = L y — 2.3L, while the perpendic- 
ular dimension was varied: (a) L z = 2.3L; (b) L z = 10.35L; 
(c) L z = 13. 8L. The corresponding free energy barriers AF 
are: (a) 1.52 ± 0.05; (b) 2.47 ± 0.13; (c) 2.29 ± 0.15, in units 
of ksT. The error bars indicate the magnitude of the scatter 
in AF for a number of independent measurements. 



during one such biased simulation. Indeed, we observe 
frequent transitions between the isotropic (52 ~ 0) and 
the nematic phase (52 ~ 1). Also shown in Fig. El is the 
corresponding time series of the reduced density. In case 
a perfect estimate for P(N) could be provided, the mea- 
sured distribution in the biased simulation will become 
flat in the limit of long simulation time. The deviation 
from a flat distribution can be used to estimate the error 
in P(N), or alternatively, to construct a better estimate 
for P(N). The latter approach was in fact adopted by 
us. First, successive umbrella sampling is used to ob- 
tain an initial estimate for P(N). This estimate is then 
used as input for a number of biased simulations using 
the modified Hamiltonian, and improved iteratively each 
time. 

To obtain the interfacial tension, the most straight- 
forward approach is to fix the lateral box dimensions 
at L x = L y , and to increase the elongated dimension 
L z 3> L x until a flat region between the peaks in the dis- 
tribution P(N) appears. For soft spherocylinders of elon- 
gation L/D = 10, the results of this procedure are shown 
in Fig. El Indeed, we observe that the region between the 
peaks becomes flatter as the elongation of the simulation 
box is increased. Unfortunately, even for the largest sys- 
tem that we could handle, the region between the peaks 
still displays some curvature. In other words, the inter- 
faces are still interacting, indicating that even more ex- 
treme box elongations are required. Ignoring this effect, 
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FIG. 7: Finite size extrapolation of the IN interfacial tension 
of soft spherocylinders with e = 2 and rod elongation L/D = 
10 and 15. Shown is the interfacial tension of the finite system 
7(L X ) in units of ksT/D 2 , measured in a cubic system with 
edge L x , as a function of (L/L x ) 2 . Lines are linear fits to 
the data using Eq.lJlJ with 6 = 0. The upper (lower) arrow 
indicates the estimate of 7in obtained using the method of 
Fig. El for L/D = 15 (10). 



and applying Eq. © to the largest system of Fig. El we 
obtain for the interfacial tension tin = 0.0022 k^T/D 2 . 
For rod elongation L/D = 15, the distribution of the 
largest system that we could handle is shown in Fig. 
The height of the barrier reads AF = 10.6 /cbT, and the 
corresponding interfacial tension 71N = 0.0053 k^T/D 2 . 

An alternative method to obtain the interfacial ten- 
sion is to perform a finite size scaling analysis. Following 
Ref. I28L the interfacial tension j(L x ) in a cubic system 
with edge L x , shows a systematic L x dependence that 
can be written as: 

l(L x ) = 700 + a/Ll + b ln(L x )/L x , (4) 

with the interfacial tension in the thermodynamic 
limit (assuming periodic boundary conditions and dimen- 
sionality d = 3). In general, the constants a and 6 are 
not known. However, recent theoretical arguments |35| 
suggest that in three dimensions, the logarithmic term 
should vanish, implying 6 = 0. To estimate 700, we 
used Eq.© to measure 7(L X ) for a number of differ- 
ent system sizes. We then used Eq.Q to extrapolate 
these measurements to the thermodynamic limit, assum- 
ing 6 = 0. For soft spherocylinders, the results of this 
procedure are summarized in Fig.0 Shown is the interfa- 
cial tension of the finite system as a function of (L/L x ) 2 . 
The data seem reasonably well described by Eq. @ , as is 
indicated by the fits. The corresponding estimates for 
the interfacial tension are 71N = 0.0035 k^T/D 2 and 
7in = 0.0059 k B T/D 2 , for L/D = 10 and 15, respec- 
tively. 
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TABLE I: Bulk properties of the coexisting isotropic and nematic phase in a system of soft spherocylinders interacting via 
Eq.Q with e = 2 and rod elongation L/D = 10 and 15. Listed are the reduced density p* and the normalized number density 
pLD 2 of the isotropic and the nematic phase. Also listed is the interfacial tension ttn of the IN interface, obtained using finite 
size scaling, expressed in various units to facilitate the comparison to other work. The error bar in the latter quantity indicates 
the uncertainty of the fit in Fig. 



L/D 


isotropic phase 


nematic phase 


interfacial tension 7in 






p* pLD 2 


p* pLD 2 


~W LD 


k B T 
(L+D)D 


10 


0.363 0.388 


0.397 0.424 


0.0035 ± 0.0003 0.035 


0.039 


15 


0.244 0.267 


0.280 0.307 


0.0059 ± 0.0001 0.089 


0.094 



For comparison, the arrows in Fig. [7] mark the inter- 
facial tension as obtained using the previous method of 
Fig. Clearly, there is some discrepancy. The prob- 
lem related to the first method is that the system size 
was not sufficient to completely suppress interface inter- 
actions. Moreover, the lateral L x dimension was also 
rather small, so there may still be finite size effects in 
this dimension. Hence, we believe the finite size scal- 
ing results to be more reliable. The latter estimates are 
listed in Table U together with the coexisting phase den- 
sities, which effectively summarizes the main results of 
this work. To our knowledge, this is the first study to 
report a systematic finite size scaling analysis of the IN 
interfacial tension. The results of Fig.[7|seem reasonable, 
but simulations of larger systems are clearly needed, in 
order to confirm the validity of Eq. Q in systems of elon- 
gated particles. The advantage of the present simulation 
approach is that the statistical errors are small, and that 
finite size effects are clearly visible as a result. 



VI. DISCUSSION 

In this section, we compare our findings to other work. 
More precisely, we consider (1) theoretical treatments 
within the Onsager approximation, (2) theoretical treat- 
ments beyond the Onsager approximation, and (3) other 
simulations. For reasons outlined in the introduction, we 
do not compare to experimental data. 

It is clear from the phase diagram of Fig. 0] that the 
Onsager limit is not recovered until for very large rod 
elongation, exceeding at least L/D = 40. As a result, 
our estimates for the interfacial tension differ profoundly 
from Onsager predictions. Typically, 71N in our simula- 
tions is four times lower compared to Onsager estimates. 
Note that our simulations also show that 71N increases 
with L/D, towards the Onsager result, so there seems 
to be qualitative agreement. However, to properly access 
the Onsager regime, additional simulations for large elon- 
gation L/D are required. Unfortunately, as indicated by 
the scatter in the data of Fig.^] and also in Ref.lT^. such 
simulations are tremendously complicated. It is ques- 
tionable if present simulation techniques are sufficiently 
powerful to extract 7n\r with any meaningful accuracy in 
the Onsager regime. 



If we compare to the theory of Ref. IToL which goes 
beyond the Onsager approximation and should there- 
fore be more accurate for shorter rods, we observe bet- 
ter agreement. For L/D = 10, the theory predicts 
7in = 0.0877 k B T/(L + D)D, which still differs from our 
result by a factor of approximately 2. For L/D = 15, 
however, a naive interpolation of the data in Ref. flOl 
yields tin ~ 0.1 k B T/(L + D)D, which exceeds our result 
by only 6%. Note that Ref. 10 considers hard sphero- 
cylinders, whereas our work is based on soft spherocylin- 
ders. The simulations of Ref. [l8| on ellipsoids suggest that 
the interfacial tension increases, when switching from a 
hard to a soft potential. The good agreement we observe 
with Ref. Hol should therefore be treated with some care. 

As mentioned in the introduction, computer simula- 
tions of soft ellipsoids with k = 15 yield interfacial 
tensions of 71N = 0011 ± °- 004 k B T/B 2 and 7 i N = 
0.016 ±0.002 k B T/B 2 HHH- For L/D = 15, our result 
for soft spherocylinders is considerably lower. Obviously, 
spherocylinders are not ellipsoids, and this may well be 
the source of the discrepancy. Note also that the shape of 
the potential used by us is different from that of Refs.ll8l 
and|U 

In summary, we have performed grand canonical Mon- 
te Carlo simulations of the IN transition in a system of 
soft spherocylinders. By measuring the grand canoni- 
cal order parameter distribution, the coexistence densi- 
ties as well as the interfacial tension were obtained. In 
agreement with theoretical expectations and other sim- 
ulations, ultra-low values for the interfacial tension 71N 
are found. Our results confirm that for short rods, the in- 
terfacial tension, as well as the coexistence densities, are 
considerably lower than the Onsager predictions. This 
demonstrates the need for improved theory to describe 
the limit of shorter rods, which is required if the connec- 
tion to experiments is ever to be made. In the future, we 
hope to extend our simulation method to the case of hard 
spherocylinders. Note that grand canonical simulations 
of hard particles are challenging, because the acceptance 
rate for insertion is typically very low. We are currently 
investigating different biased sampling techniques in or- 
der to improve efficiency. Also the investigation of the 
structural properties of the IN interface is in progress. 
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